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Preservation of Takens-Bogdanov bifurcations for delay 
differential equations by Euler discretization* 

YiNGXiANG Xu^and Chengchun Gong* 



Abstract 



A new technique for calculating the normal forms associated with the map restricted 
to the center manifold of a class of parameterized maps near the fixed point is given first. 
Then we show the Takens-Bogdanov point of delay differential equations is inherited by 
, the forward Euler method without any shift and turns into a 1:1 resonance point. The 

ySL • normal form near the 1:1 resonance point for the numerical discretization is calculated 

I next by applying the new technique to the map defined by the forward Euler method. The 

(~| . local dynamical behaviors are analyzed in detail through the normal form. It shows the 

Hopf point branch and the homoclinic branch emanating from the Takens-Bogdanov point 
are 0(e) shifted by the forward Euler method, where e is step size. At last, a numerical 
experiment is carried to show the results. 



^ ■ 1 Introduction 

^T) . Numerical methods may take on many possible phenomena when applied to certain dynamical 

CO ■ systems. It is of great importance to investigate what kind of properties of the original systems 

~~~J^ , could be preserved by discretization. Tremendous researches on numerical stabilities of sorts 

CSJ ■ are the first investigation of this type where various sufficient and necessary conditions are 

, developed for numerical schemes reproducing the asymptotic stability of differential equations, 

■ one can refer to the monographs [5j and [2, for ordinary differential equations (ODEs) and delay 

differential equations (DDEs), respectively. Whether or not the numerical discretization will 
inherit the bifurcations of the original systems is another important research field of this type. 
, In this paper, we consider the DDEs of the type 

z{t) = f{z{t),z(t-l),a), (1) 

where z G K", a S is a bifurcation parameter, /(zi, Z2, a) is C""(r > 2) smooth with respect 
to zi,Z2 and a. The state space of ([U, denoted by C = C([— 1, 0], E"), is a Banach space of 
continuous mappings from [—1,0] to M" equipped with norm ||(/)|| = max0g[_i q] 10(^)1 (I ' I is 
some norm in M"). The equation ([Ij is assumed to undergo a Takens-Bogdanov bifurcation 
near {z, a) = (0, 0). Then, in the parameter plane (ai, there exist a Hopf point branch and 
a homoclinic branch emanating from the Takens-Bogdanov point of equation see [35]. In 
this paper we show that the forward Euler discretization, when applied to equation ([IJ, can 
preserve the bifurcation structure near the Takens-Bogdanov point (z, a) = (0, 0) of ([T]) by an 
0{e) shift, where e is the step size of the Euler method. 
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Problems of this type have been focused on ODEs for many years. The Hopf bifurcation 
accepted the most attentions in this field after the work of Hofbauer and looss |15j, where 
they proved the forward Euler discretization exhibits the Hopf bifurcation of the same type as 
the continuous system undergoes. In 1998, Wang, Blum and Li ^19j gave the most complete 
results for the codimension 1 bifurcations: the general one step methods of order p and specific 
methods like Euler, backward Euler, explicit and implicit Runge-Kutta methods are proved to 
inherit the elementary bifurcations, such as saddle-node, fold, pitchfork, Hopf and transcritical 
bifurcations, of the continuous systems. For connecting orbits, Beyn [3] first showed the exis- 
tence of homoclinic orbits under discretization. Fiedler and Scheurle [131 , Zou and Beyn [25] 
proved the general one step method exhibits the discrete connecting orbits approximating to the 
one of the original ODEs by the order of the method, independently. Most recently, the major 
concern for this problem is whether the numerical discretization can inherit the codimension 2 
bifurcations of ODEs. Loczi and Chavez [18] showed the Runge-Kutta method will reproduce 
the nondegenerate fold, cusp and Takens-Bogdanov singularities of general ODEs without any 
shift. For generalized Hopf bifurcations, the Hopf point branch emanating from the generalized 
Hopf point is shifted by the order of the general one step method used as expected, however, 
the generalized Hopf point is shifted of the first order, regardless the order of the numerical 
method used, and turns into generalized Neimark-Sacker points. This result seems not nat- 
ural, but the numerical experiments shows a much better result cannot be expected, cf. [7]. 
For the Hopf-Hopf bifurcation, the one step method will reproduce the bifurcation point as a 
double Neimark-Sacker point by an 0{e^) shift, as well as the Hopf point branch bifurcated 
from the point, see Chavez [6], under the assumption that the ODEs undergo the fold- Hopf 
or Takens-Bogdanov bifurcation, proved the singularities will be persisted by the Runge-Kutta 
method and the Hopf point branch emanating from the singular points is 0{eP) shifted (e is step 
size). In addition, the discrete fold point branch and the discrete Hopf point branch emanating 
from the singularities is shown to intersect transversally. However, according to what we have 
seen up to now, no discussion about the preservation of the homoclinic branch emanating from 
the Takens-Bogdanov point by the numerical scheme is available. 

In these days, the preservation of dynamical behaviors for DDEs by numerical discretizations 
received many considerations as well. Wulf and Ford first showed the forward Euler method 
will reproduce the same type Hopf bifurcations as the scalar DDE with one constant delay 
undergoes [21| . The subsequent numerous relative works had led this problem to a more general 
case of both the numerical methods (Runge-Kutta method, linear multi-step method etc.) and 
the type of the differential equations (with multiple delays or state-dependent delay, etc.). 
For other dynamics of DDEs preserved by discretization, we recall the following results. Liu, 
Gao and Yang [17] proved the Runge-Kutta method preserves the oscillations of the equation 
x{t) + ax{t) + aix{[t — 1]) ~ 0. The asymptotically stable periodic orbits of the autonomous 
DDEs with one constant delay were proved by In't Hout and Lubich [16] to be shifted by 
the s-stage Runge-Kutta methods with the discretization order p. For a more general delay 
equation 

i = f{xt), (2) 

where Xt{ff) = x{t + 6), Farkas proved by explicit calculation that the unstable manifold will 
be close to its Euler discretized counterpart if the step size e is sufficiently small [9 as well 
as a numerical shadowing result which obviously means the stable manifolds of ([2]) could be 
preserved numerically [HI. Xu and Zou [53] extended the results of [IS] to DDEs, they showed 
the homoclinic orbits should be preserved by the forward Euler method with a shift of 0{e). 
In addition, the possibility of extending the numerical scheme to the implicit method and the 
Runge-Kutta method is discussed there. 

The problem we care most is whether the bifurcation diagram near Takens-Bogdanov point 
of DDEs could be inherited by numerical discretiztation. This requires us to investigate the 
dynamical behaviors of the map defined by the numerical scheme. The normal form analysis is a 
useful tool for accomplishing this aim. It allows us to know how the normal form coefficients and 
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the generalized eigenspace of the numerical scheme are related to their continuous counterparts. 
Consequently, the dynamics of the numerical scheme could be studied in detail in terms of the 
coefficients of the original equation and then we know how the bifurcation diagram is preserved 
by direct comparison. Trying to fulfill our goal in this way, we will also face some problems. 

Generally speaking, the numerical scheme for ODEs could be regarded as a map. For 
example, the one step method for solving 

X = fix), a; € R" 

has the general form 

Xk+i = Xk + e<p{xk,s), (3) 

where e is step size. Obviously, equation is a map from K" to R" and the method of 
normal forms for finite-dimensional ODEs could be modified to analyze the local bifurcation 
behavior. Where the representations of the complete eigenvectors are required to construct a 
transformation which lead the map ([3]) to a standard form could be dealt with easily (see |20] 
section I.IC for detail). In addition, the technique requires computing the center manifold first 
before evaluating normal forms for the map on the center manifold. However, the forward Euler 
method for solving 

x{t) = ax(t) + bx{t - 1), X e M" 

has the form 

Xk+i = Xk + sauk + ebxk-m, (4) 

where e = is step size, to is a positive integer. The map defined by Q is, obviously, from 
^{m+i)xn j|j(m+i)xn g^^^ -j-j^g dimensionality tends to infinity as the step size tends to 0. As 
a result, the amount of the eigenvectors is very large as the step size is small, and the method 
of normal forms for maps described in [20^ is not efficient any more. 

Highlighted by the normal form analysis for retarded functional differential equations pro- 
duced by Faria [TlJ [12] , a new method for calculating the normal forms associated with the 
map restricted to the center manifold tangent to the invariant manifold of the lienarization of 
the map of the type mentioned above at the fixed point is given. It has many advantages: one, 
without computing beforehand the center manifold of the singularity; two, allows obtaining the 
coefficients in the normal form explicitly in terms of the considered map, therefore, one can 
obtain the accurate dynamics through analyzing the obtained normal form, and three, no need 
of computing the eigenvectors associated with the noncritical eigenvalues, and then the amount 
of calculations are reduced remarkably. 

The rest of the paper is arranged as follows. We first go deep into the techniques for normal 
form calculations for the maps similar to (|4]) in section 2. The Takens-Bogdanov bifurcations 
for equation ^ is recalled in section 3. In Section 4, we first show the Takens-Bogdanov point 
of DDEs is inherited by the forward Euler method exactly and turns into a 1:1 resonance point. 
Then the normal form near the 1:1 resonance point for the numerical method is calculated by 
applying the techniques developed in Section 2. The local dynamical behaviors are analyzed in 
detail through the obtained normal form. It shows that the Hopf point branch and the homo- 
clinic branch bifurcated from the Takens-Bogdanov point of ([ij are inherited by the forward 
Euler scheme by a shift of 0{e). A numerical experiment for a scalar DDE is carried to show 
our theoretical results in the last section. 

2 Normal forms for a class of maps with parameters 

We consider a class of maps of the following type 

u 1-^ C{a,e)u + F{u,a,e), (5) 
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where u e M("+i)"' with n,m,l E e = — is a parameter used in later apphcations where e 
corresponds to the step size, a = (ai, • • • , Up) e is a bifurcation parameter. F is C^{r > 2) 
smooth with respect to u and a, and satisfies F(0, 0, e) = and §^(0, 0, e) = 0. From the map 
(O, we see that the dimensionahty of ([S]) wiU tend to infinity as the parameter e tends to 0. 

Denote C(e) = C(0, e) and F{u, a, e) — C{a, e) — C(0, e) + F{u, a, e), then we can write the 
map ([5]) as 

ui^ C{e)u + F{u,a,e). (6) 

We drop the bar from F if no confusion occurs. For simphcity, we wiU drop the dependence on 
e in the rest of this section. 

We reformulate the map ([S|) by regarding a as a variable to the following map without 
parameters 

"^.^cH + f^^"'")), (7) 



a J \aj \ Op / 
where C — diag{C, /p}. To linearize the map ^ at {u,a) — (0,0) we obtain 

which has the characteristic equation given by 

dct(A/ - (7) = ^ (A - 1)P det(A/ - C) = 0. (9) 

As seen in the former equation, we often omit the subscript of the identity matrix / if no 
confusion occurs throughout this paper. 

Denoted by A = {A e a-(C) : |A| 1} and A = AU {I,-- - ,1}, 0j the (generalized) 

p times 

eigenvectors of C associated with eigenvalues A^, i = 1, 2, • • • , (m + l)nl, then P — span{0i : 
Xi E A} is the invariant space of C associated with A. If c is the number of the eigenvalues of C 
in A, counting multiplicities, then we have dim P — c. If we denote by $c = {(ki i ' ' ' i '/'c) a basis 
for P, = (V'f , • • • , ^cV a basis for the dual space P* in m(™+i)"'*, then (^-c, $c) = h, the 
identity matrix in M.'^, where the dual form takes the scalar product of vectors. We denote by J 
the c X c constant matrix such that C^c = ^c^; its spectrum coincides with A. Denoted by Q 
the invariant space of C associated with a{C) \ A, then Q = span{(/)i, i = c+ 1, • • • , (m + l)nl}. 
If we denote $su = ('/'c+i, • • ■ , 0(m+i)ni) a basis for Q, "ifsu = ("^f+i, • • • , i^fm+i}niV a basis for 
the dual space Q* in j^l'^+i)"'*^ then we have {"^suj^su) = I{m+i)ni-c, the identity matrix in 
^((m+i)ni-c)xam+i)ni-c) ^ dcnotc by Cq the constant matrix such that C$™ = ^suCq; its 
spectrum coincides with a{C) \ A. 

With the notations above, the generalized eigenspace P of C associated with A satisfies 
P — span<i>c, where 4>c — diag($c,^p)- The left generalized eigenspace of C associated with A 
is spanned by = diag(^'c, Ip). Denoted by J the (c+p) x (c+p) matrix subject to C^c = ^cJ, 
we have J = diag( J, Ip) and its spectrum coincides with A. Clearly, cr(C) \ A = cr{C) \ A. If we 
denote Q the complementary space of P in the space k(™+i)"'+p^ then Q is spanned by ^su, the 
generalized eigenspace of C associated with cr{C) \ A. Obviously, ^su = i^'^u^^p)'^ ■ We denote 
by ^'su the left generalized eigenspace of C associated with o'(C) \ A, then ^'s^ = {'^su,Op)- 
Obviously, we have (*s„,l>s„) = I(,n+i)ni-c- 

According to the decomposition of the state space R(™+i)"' = PQ)Q, u could be represented 
as M = <^cx + <^suy with X EW,yE 

Therefore, projecting map ([7|) onto P and Q respectively, we get 



a J \aj \ Op J' (10) 

CQy + ^suF{^cX + <S> suU-i ■ 
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Noting that a is mapped into itself, the previous map exactly is 



where a should be considered as a variable. In addition we recall that the matrices J and Cq 
are of the form of uptriangular. 

The first part in map (jlip has a fixed dimensionality, that is c, the dimensionality of the 
center manifold of ^ associated with A. While the second part has a variable dimensionality, 
which tends to infinity as the parameter e tends to 0. This makes it very hard to give explicitly 
the expressions for the (generalized) eigenvectors of ([8]) associated with (j{C) \ A. Besides, 
the amount of the eigenvectors will also tend to infinity as the parameter e tends to 0. It is 
impracticable to get the map restricted to the center manifold of ([TT|l by the technique described 
in |20j directly. In the next we go deep into develop a method which allows getting the normal 
forms associated with A on the center manifold avoid computing the eigenvectors '^su and ^ su- 

Taylor expanding -F(u, a) with respect to {u, a) leads to 



F{u,a) = ^iF,(u,a), {u,a) G r(™+iW+p. (12) 
Defining/, = (/I, /2) with 



J- 

j>2 



fj{x,y,a) = *c^j($ca; + $suy,a), 
ff{x,y,a) ='^suFji'S>cX + 'S>suy,a), 



(13) 



the map ([T]) is equivalent to 



x^ Jx+ -^f}ix,y,a), 
y-^ CQy+ E —ff{x,y,a) 



(14) 



>2 J 



where xeR" andy e rC^+i)"'-^ 

The normal forms are obtained by a recursive procedure, computing at each step the terms 
of order j > 2 in the normal form from the terms of the same order in the original map and 
the terms of lower orders already computed for the normal form in previous steps, through a 
transformation of variables ^ 

{x,y) = {x,y) + —Uj{x,a), (15) 

with X, xe y, ye and = (Uj, Uf) G t//+P(M'=) x i/'=+P(]R(™+i)»i-c), where, 

for a normed space X, Vj^'^^X) denotes the linear space of homogeneous polynomials of degree 
j in c + p real variables, {x,a) — {xi,X2, ■ ■ • , Xc, ai, • • • ,a.p), and with coefficients in X, 

V^^+^iX) = { C(,,)x''a' : (g, /) G Ng+^ C(,,) € X}, 

l(9,i)l=J 

x'^a' = xl^ ■ ■ ■ xf'a-l ■ ■ ■ dp for g = (gi, • • • , qc) G Wq and I — {li, - ■ ■ ,lp) e Nq, with the norm 



l(9,01=i 



I] |C(9,0U- 

\{q,i)\=o 
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We assume that after computing the normal form up to terms of order j — 1 the map is 



Jx+ J2 -.gli^iV^OL) + —J}{x,y,a) 
i=2 j! 



1 2 1 

where 



y ^ CQy+ E -^9tix,y,a) + —ff{x,y,a) + ■ 



(16) 



= fjix,y,a) - [C/j(Ja;,a) - j;7/(x,a)], 
g]{x,y,a) = ff{x,y,a) - [Uf{Jx,a) - CqUfix^a)]. 
These formulas can be written for gj — ( g j ,g'j) as 

ff, = /,-M,C/„ (17) 

with fj = ifjifj) ^^^d Afj- defined below. 

Definition 2.1 For j > 2, let Mj denote the operator defined in y/+P(]R'= x ^('"+1'"'-'=), with 
values in the same space, by 

M,{p,h) ^[M]p,M]h) 
{Mjp){x,a) = p{Jx,a) — Jp{x,a) 
{Mjh){x,a) = h{Jx,a) ~ CQh{x,a), 

with domain D{M.j) = y/+P(M^) x 

Then the problem we met is how to choose Uj so that gj has a simple form. 
We decompose the spaces y/+^(M^) and v^'=+p(m(™+i)»'-c) as 

V-^^\W) = Im(Mj) ® Im(M/)^ 
V^+P{W) = keriMj) ® ker(Mj)=, 

and 



. ^ ^^^(rn+i)«i-c) ^ lm{Mf) ® Im(M|)S 
(«("'+!)"'-'=) = ker(M2) © ker(M2)^ 

We denote the projections associated with the above decompositions of V^'^^^{M.'^) x ]/'^+''(m('™+i)"'^^) 
over Im(A/j^) x Im(M|) and over ker(Mj^)'^ x ker(M|)^ by, respectively, P/.j — i^ij^^ij) 
and Pkj — (PkjTPKj)- The complementary spaces in the above decompositions Im(A'fj)'^, 
ker(M*)'' (i = 1,2), are not uniquely determined. As a consequence, normal forms are not 
unique, and depend on the choices of Im(Mj)'^ (z = 1, 2). 

Let us now consider the right inverse of Mj with range defined by the spaces complementary 
to the kernels of Mj {i = 1, 2), namely M'^ = {{Mj)-\{Mf)^^) with Mr^oPijoMj = Pkj- 

Taking y = in formula ^T7\ . an adequate choice of Uj, 

Uj{x, a) = M-^Pi,J,{x, 0, a), (18) 

allows taking away from fj its component in the range of Mj , leading to 

g,{x,G,a) = {I-Pi,,)f,{x,0,a). (19) 

Therefore, the normal form for map ([5]) relative to the invariant space P and the projections 
Pi J, Pk,j (j = 2, 3, • • • ) is the map in K'^ x 

X ^ Jx + J2i>2 -J9j{x,y,a) 
y ^ Cqy + Ej>2 —^9]{x,y,a). 
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In the next, we show that the normal form relative to P takes a more simple form of 
c-dimensional map on a locally invariant manifold for ([5]) tangent to P at zero. 

Recall that the matrices J and Cq are of Jordan forms, then by the Theorem 3.8 in [1] we 
know that the normal form (|20|) can be chosen so that its nonlinear part contains only resonant 
monomials. 

Definition 2.2 We say that the map ^ satisfies the nonresonance conditions relative to A e 
ct{C), if 

A« ^ fi, for all ^ e a{C) \ A, q G Nq+p, \q\ > 2, 

where A — (Ai, • • • , Xc+p) and Ai, • • • , Xc+p are the elements of A, each one of them appearing 
as many times as its multiplicity as an eigenvalue of the matrix C. 

Clearly, the fact that the map ([7]) satisfies the nonresonance conditions relative to A e cr(C') 
is equivalent to the map ([5]) satisfies the nonresonance condtions relative to A g o'(C) since the 
difference between A and A are the eigenvalues introduced by regarding a as a variable, that is 
A\A = {V_^}. 

p times 

Theorem 2.3 Consider the map and let P be the invariant subspace ofC associated with 
a nonempty finite set A of eigenvalues. For the decomposition o/M'™^-'^)"' = R^©R(™+-'^)"'~'^, we 
get u = ^cX+^suU, where x £ M"^ and y £ ^{m+i)ni-c ^ formed by the generalized 

eigenvectors of C associated with A and cr(C) \ A, respectively. If the nonresonance conditions 
relative to A are satisfied, then there exists a formal change of variables {x,y) i— >■ {x^y) of the 
form X = X + p{x,a),y = y + h{x,a), such that the map ^ is equivalent to the following map 
m W X R(™+i)"'-c 

X ^ Jx + J2j>2 —,9]{x,y,a) 

- (21) 
y ^ Cqy + Y.j>2 j^9o{x,y,a), 

where g],g^j are computed as with g'^(x,0,a) — for all j > 2. This map is in normal 

form relative to P. If there exists a locally invariant manifold for the map tangent to P at 
zero, then it satisfies y ~ and the map on it is given by the c-dimensional map 

x^ Jx + J2^g^{x,0,a), (22) 
which is in normal form for maps. 



3 Takens-Bogdanov bifurcations of parameterized delay 
differential equations 

To show the bifurcation structure near Takens-Bogdanov point could be preserved by numerical 
discretization, we recall the results on the bifurcation structure near Takens-Bogdanov point in 
this section, cf. [22]. 

We consider the DDEs of the type 

z{t) = fiz{t),z{t-l),a), (23) 

where z £ R", a € is a bifurcation parameter, f{zi,Z2,ct) is a C^{r > 2) smooth function 
from R" X R" x ^ ^^^^ 

/(0,0,a) =0,1^(0,0, a) = 0,1^(0,0, a) =0, Va G M". (24) 

OZi OZ2 
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Denote A = ^-(0, 0, 0) and B = ^■(0, 0, 0), then we can rewrite equation ((23)) as 

z{t) = Az{t) + Bz{t - 1) + F{z{t), z{t - 1), a), (25) 
which could be hnearized at (z, a) = (0, 0) as 

z{t)^Az{t)+Bz{t-l). (26) 
The characteristic equation of (1^51) is given by 

det(^/ - A - Be"^) = 0. (27) 
We assume {z,a) = (0,0) is a Takens-Bogdanov point of (|23|) . i.e., 

A is a root of (P7)) with algebraic multiplicity 2 and geometric multiplicity 1, and the other 
roots exhibit nonzero real parts. 

Lemma 3.1 [22] There exist e M"\{0}, 0^ E W\ and ipl e M"*\{0}, Vi e K"*; ^^ej/ satis/j/ 
the following equations 

(1) (A + B)^^ = 0, (2) {A + B)^l = {B + 1)^1 

(3) V2°(^ + i?) = 0, (4X(^ + i?) = V^0(i? + /), 

(5) ^-2^^ - h^'2Bcb"i + ^^Bd^l = 1, (28) 

(6) - + ^^^B4>1 + i^'2°S0? - \^IB4>1 = 0. 

Taylor expanding F(z{t), z{t — 1), a) with respect to z{t), z(t — 1) and a we obtain 

F(z(<), z(t - 1), a) = ^ ^F,{z{t),z{t ~ 1), a), (29) 

j>2 

where the first term (j = 2) can be expressed in the form 
^F2iz{t),zit-l),a) 

= Aiaiz{t) + A2a2z{t) + Biaiz{t - 1) + B2a2z{t - 1) (30) 
+ Er=i EMt)^{t - 1) + Eti F^z,it)z{t) 

+ lT^=lG^Z,{t~l)z{t~l) 

with Ai,Bi{i — 1,2), Ei, Fi,Gi{i = 1,2,--- ,n) coefficient matrices, and there is no terms of 
0(a2) in F2{z{t),z{t - l),a) since F(0,0,a) = 0,Va € K^^ 
Denote 

a= ^2° Er=i(^«+^«+G.)0? 



and 



with 



'^>n"M (32) 

K2/ V"2 ' 



n 



V'2°(Al+Sl)0? V2°(^2+S2)</'? \ 

{<(Ai+i3i)</.? {^1{A2+B.2)cj>l 

+V^O((^l+Sl)0§-Sl0?)} +^2°((^2+S2)0§-i?2</'?)} y 

For the bifurcation structure near the Takens-Bogdanov point of DDE we have the fol- 

lowing theorem. 
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Theorem 3.2 \22\Assume the assumption A holds, detll ^ and a ■ b =/= 0. Then there 
exists a constant > 0, such that when < Ki{ai,a2) < Ki, in the parameter plane (cki,ck2) 
near the origin there exist two curves: Ih and loo 

1. the curve Ih, which has the following local representation: 

lh = {{ai,a2): K2(ai, 02) - -Ki(ai, 02) + /i.oi. = 0, Ki(ai, a2) > 0}, 

a 

is a H op f point branch oj the DDE where h.o.t. = o(|(ai, 02)!), i.e. Ih consists of Hop f 

bifurcation points of [^) : 

2. the curve l^o, which has the following local representation: 



^00 = 02) : h{ai,a2) + h.o.t. = 0, £12) > 0}, 



(33) 



is a homoclinic branch of the DDE (23\), where (01,02) = '«2(cki, Ck2) ~ >^i{ctii <^2)) 
Ki (01,02), /i(-) is a continuously differentiable function with /i(0) — |6a^^ and h.o.t. — 
0(1(01,02)!). In other words, equation h23\) has a unique homoclinic orbit connecting the origin 
for each (01,02) G loo- 



4 Preservation of Takens-Bogdanov bifurcations by Euler 
discretization 



, m G Z+ is given by 



The forward Euler scheme for solving (|23p with step size e = 

Zk+l = Zk+ ef{zk, Zfe-m, o), 

which can be reformulated to 

Zfe+i = Zfe + eAzk + eBzk-m + £F{zk, Zk-m,a). 
Denoting Uk — {z]: , z]^_^, . . . , Zk_jn)'^ G ([SS]) can be rewritten as 

Uk+i = Cuk + H{uk,a), 



(34) 
(35) 
(36) 



where 



C 



( I + eA 
I 





V 



Q eB \ 








and 



0)^,0, 



/ 



(37) 



H{uk,a) = {eF{zk,Zk 
Linearizing the map (pS)) at (u, o) = (0, 0) we obtain 

Uk+i = Cuk, 
which has the characteristic equation 

det(A/ - C) = 0. (38) 

Noting (|36p is a map of the type of ([5]) as Z = 1, the techniques developed in Section 2 allow 
us to get the normal forms on the low dimensional center manifold without computing the 
eigenvectors associated with the eigenvalues of the modulus other than 1. While the map ([5]) 
as Z > 1 could correspond to more general numerical schemes, eg. the general one step method. 

Under the assumption that undergoes a Takens-Bogdanov bifurcation at (z, o) — (0, 0), 
we first show the Takens-Bogdanov point (z, o) = 0, 0) of is inherited without any shift by 
the Euler method (|35|) and turns into a 1:1 resonance point. 



9 



Theorem 4.1 Assume the assumption A holds. Then ^34\ ) undergoes a 1:1 resonance at 
(z,a) = (0,0). 

Proof. We only need to show the double zero eigenvalue of (pS)) turns into a double unit 
multiplier, \i 2 = 1 at a = 0, and no other eigenvalue of p7p has modulus 1. 
Assume 

=0, (C-/)02 -01, 
MC-I) =0, MC-I) =^2, 



(39) 



and require 

We can check that the following choice of 0i, 02, "02, "01 

0i=e(0r,---,0rr, 

02 = £{m(t>S^ , m(t)^ — 0°"^, • • • , rn02"^ — mctf^Y" ^ 

^2 = i_ Aoj3w,o (V'2°,g02"i?,--- ,£^2°i?), (40) 
fulfill the requirements above. Besides, the Fredholm alternative Theorem implies 

V'2 • 01 = '01 • 02 = 0. 

These show that 1 is an eigenvalue of p7p with algebraic multiplicity 2 and geometric multi- 
plicity 1. 

Next, we show no other eigenvalue of p7p has modulus 1. 

Denote dip) = e^'{^J.I - A) - B and D{^i,e) = e^'{g{pe)^J.I - A) - B, where 5(2;) 
Then, see we have det d{fi) = is equivalent to (l?71) . and if det D{fi^,£) — 0, then A = e^^^ 
solves (1551) . Besides, we have limZ)(/x,e) = d(ii). Hence for any series fi^ of solutions of 

detD(^,e) = 0, there exists a solution /xq of det c?(/i) = such that lim ~ ^q. Therefore, 

if other than the eigenvalue 1 there exists any eigenvalue of denoted by A^, s.t. |Ae| = 1, 
then = i In A^ solves det -D(/i, e) = 0. Since ^.^ is pure imaginary, we have /xg, as the limit of 
/ig as e — >■ 0, is either pure imaginary or and solves (j27p . This contradicts to the assumption 
A. □ 

Let $c — (01,02), = ("01,02), and A = {1, 1}, then we know the invariant space P of 
(j36p associated with A is spanned by $c, the dual invariant space P* of (|36)) associated with A is 

1 1 



spanned by ^-c. They satisfy (*c, $c) = h, C^c = ^cJ, and *cC = J*c with = 1 g ^ 

We then consider the Taylor expansion of H{uk,a) with respect to {uk,a)- By (P^ . we 
have 

H{uk,a) ^ e{F{zk,Zk-.m,af ,0, - ■ ■ ,0f = (^ ^F,-(zfc, Zfe-™, a)^, 0, • • • ,0)^, (41) 

j>2 ^' 



where 



= AiaiZk + A2a2Zk + BiaiZk^m + B2a2Zk-m (42) 
+ ELi E,zlzk-m + Er=i F^ZkZk 
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Evidently, the canonical basis of V2'*(R^) is composed by the elements 



fxiX2\ /a;iai\ /a;iQ;2\ /a;2ai\ 
[ Q J' [ J' [ )'[ J' 



fx2a2\ faia2\ 

{ )•[ )' 



i ° 

\XiX2 





Xiai 







/ 









aia2 



\X1a2J' \X2aiJ' \X2a2^ 

and the images of these elements under are, respectively 



/ 2x1X2 + x: 

\ 







2X12:2 + Xo 



xl\ fx2ai\ fx2a2\ fO 

oj' j' )\q 

XlX2\ 



— xiai\ 



-a;ia2 

X2a2 



-a:;2ai 
/' 



\ /-X2a2\ (-axa2\ 



Therefore, a basis of lni(Af2 can be taken as the set composed by the elements 











2 / ' I „,2 ; ' 







/ \ / 



\X\X2) V^i'^i 





a;ia2 



\X20L\ 







( 



X2OL2) \0i\0i2 



Denoting by $2 the first n rows of <i>c, the last n rows of $cj we have i?($ca 
Denoting by the z-th element of , we have from (j30|) 



1 



-F2[[e^\Al)xAe'KA2~^\)x,a) 



.o^. ...0 ^0 

^iai(e(/)?,(/)^)(xi,X2)^ + A2a2(£0?,(/)2)(X1,X2)^ 

+Biai(e(/.;, (^2 - <fi){.xx,X2f + B2a2(£0?, </)^ - <\>\){xx,X2f 



+ {e4>l Alr){xi,X2f {e4>l ,4>l- 4>l){x^,X2)^ 

n 

+ Y.F,{e4>l,,<tP2i){xuX2f{e4>l4>l){xuX2f 

n 

+ G^{£^l^,^2^ ^ 0?.) (^1 , 22)^ , 0^ - 0?)(xi, 22)^ 

1=1 

e{Ai + Bi)(t>laixi + e{A2 + B2)?!)? 0:22:1 

+ + - Si0?)aiX2 + ((^2 + B2)<^^ - B24>\)a2X2 

n 

n 

+e5^{(i?, + + G,)(0?,(/.° + (/.°,(/.;) - (i?. + 2G,)4>Ul}xiX2 

n 
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Base on the expansion above, and the canonical basis of V2'(K^), Im(Af2) and Im(M2 )'^, not- 
ingthat/i(a;,0,a) = ^',i72($ca;,a) - ^(ir^n^p^V'?, TT;n|p^V'2°)^2(M, <^^)a:, (e^?,^^ 

0]')x,a), we can compute the function g2{x,0,a) = (/ — P/2)/2 (^^j "^)- By Theorem I2.3| we 
have the following results. 

Theorem 4.2 Assume the requirements in Theorem \ are fulfilled. Then the numerical 
scheme ^34^ will undergo a 1:1 resonance at (z, a) = (0,0). In addition, the numerical scheme 
fff^l) could be reduced to a 2 dimensional map on the center manifold at {z,a) — (0, 0) as follows 



Xi 1-^ Xi + X2, 

X2 ^-^ X2 + k\xi + K2X2 + a^x\ + ¥xiX2 + h.o.t., 



2,^e._.. (43) 



where ^ 



wit/i Ki,K2 a'l'^ 0,^ defined in iS^) and iSl]) . respectively. 

It is known that for the reduced map (|43l) . if • 6^ 7^ (equivalently a ■ b ^ 0), the local 
bifurcation structure near (x, a) = (0, 0) is determined by the linear and quadratic terms, and 
not the terms of order higher. Hence we turn to investigate the local bifurcation structures of 
the map 

Xi l-> Xi + X2, 

X2 '-^ X2 + nfxi + K2X2 + a'^xf + ¥xiX2 



(44) 



Lemma 4.3 Let A^(Kf , kI) be the eigenvalues of the Jacobian of |^^| ) at (— 0). Then, whe 
■¥ ^0 andO < kI <2 



\XfK,4)\^\ll + n-2 

¥ 

Hence we conclude that each point on the line segment /| = {('*ii'*2) • ^^2 ~ — '^i ~ '^iiO < 
Kf < 2} in the parameter plane (k^jKj) is a Neimark-S acker bifurcation (also known as the 
Hopf bifurcation for map) point of the map \44^ . 

Proof. Evidently the map has two fixed points, (—-^,0) and (0,0). Direct computations 
show the eigenvalues of the Jacabian at {—-^, 0) reads 



A^('^!, «5) = ^(2 + «:"2 - -,^\) ± V (2 + 4 - - 4(1 + 4 - + «f )> 

2 V a*^ 

which have a modulus given by 



Therefore, when {k\, k!^) changes from one side of /| to the other, in the parameter plane the 

¥ 

eigenvalues X^{n\, will cross the unit circle from outside to inside (kj + > 0) or 

from inside to outside {k\ j'*! + '^i < 0). Let 9 — arctan ^ ^^^^ ^ — , obviously we have 

when < Kf < 2 

e^^^V 1, for A; = 1,2,3,4. 
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Applying of the Neimark-Sacker bifurcation Theorem implies X2, , k|) — ( -,0,kI, — k!- 

T-r-TT ^'^ '^^ 

) is a Neimark-Sacker bifurcation point of the map as < < 2. In fact, there is a 
neighborhood of ( — ^,0) in which a unique closed invariant curve bifurcates from ( ^ q) 

when k| ^Kf + Kf changes signs. In other words, the line segment Z| in the parameter plane 

(^f , kI) is a Neimark-Sacker point branch of the map (|44|) . □ 

Remark 4.4 Noting the expression of w e^Ki, t/ie statement of < 2 requires the step size 
e of the numerical scheme \3J^ should be taken nicely small to reproduce the Hopf bifurcations 

ofm- 



In the next we consider the homoclinic curves bifurcates from the fixed point (x, a) = (0, 0). 
Applying the transformation of Xi = a'^{xi + j^), X2 = 0^X2 to (|44|) leads to another typical 
normal form for 1:1 resonance, that is 

Xil->-Xi+X2, 

X2^ X2~ + (kI -~ -2^)X2 +Xi + -^XiX2- 

Applying the results in [4 to the former map shows that the homoclinic curves bifurcated from 
the fixed points (x, a) = (0, 0) is depicted by 

nl^~n\-^-K\ + 0{{K\)i) (46) 

as < k\ < k", where is some positive constant. 

Remark 4.5 In fact, there exist two curves q;J(q;i) and 012^ (ai) respectively corresponding to 
the first and the last homoclinic tangency, they are exponentially close to one-another. If a is 
located in the region confined by these two curves, the map |^^[ ) possesses transverse homoclinic 
trajectories, see for detail. But this is not our goal in this paper. 



Noting that the map (|44)) is locally topologically equivalent near the origin to (|43l) , we have 
the following bifurcation results based on Lemma [4.31 and the discussions above. 



Theorem 4.6 Assume a"^ ■ b"^ ^ (a'^ ,b^ G K defined in Theorem \4-S^ - Then there exists 
a constant k^^ = niin{2, k^} > 0, such that when < Kf (ai, 02) < Ki^, in the parameter plane 
{ai,a2) near the origin there exist two curves: l\ and 1%^ 

1. the curve If^, which has the following local representation: 

Ifi = {{<^i,0i2) ■ K^io^iiOi^) - — Kf(ai,a2) + kI + h.o.t. = 0, < Kf(ai,a2) < k?^}, 



is a Neimark-Sacker point branch of the numerical scheme jS^^ , where h.o.t. = e ■o{\{ai,a2)\), 
i.e. Z| consists of Neimark-Sacker bifurcation points of {34^ ; 
2. the curve If^, which has the following local representation: 

llo = ■ h{ai,a2) + h.o.t. = 0, < K,i{ai, 02) < >i(^}, (47) 

is a homoclinic curve of the numerical scheme {^p, where h{ai, a2) = '^2~f ^'^i + f '^i; h.o.t. = 
■ o{\{ai, Q?2)|). In other words, the numerical scheme ^34\ ) presents a unique homoclinic curve 
connecting the origin for each (ai,a2) S lfx>- 

Comparing Thcorcm l4.6l to Theorem l3.21 incorporating Theorem 14.11 we obtain the following 
result. 
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Theorem 4.7 Assume the assumption A holds, that is the DDE I123\) exhibits a Takens- 
Bogdanov bifurcation at {z,a) — (0,0). Then the Takens-Bogdanov point of i23\) is inherited 
without any shift by the forward Euler scheme {34^ and turns into a 1:1 resonance point. 
Moreover, there exists an Eq > 0, such that as e < Eq, the forward Euler scheme will 
reproduce the Hopf point branch and the homoclinic branch of the DDE \2S\) with a shift of 
0{e) in parameter plane {ai,a2), specially we have 

\l^-ll\ = 0{e), 

\loo-ll.\^0{e). 

5 Numerical example 

In this section we present a numerical experiment to illustrate the theoretical results. 
We consider a 1-dimensional DDE as follows 

i{t) = (1 + ai)z(t) - (1 + a2)z{t - 1) + ]^z{t)z{t - 1). (48) 

It is easy to show {z^a) = (0,0) is a Takens-Bogdanov point of pS)) . cf. [TH[52]. The forward 
Euler method for solving it is given by 

Zk+i = Zk+ e{l + ai)zk - £(1 + a2)zk-m + ]^ezkZk-m (49) 

with the step size e = — , m G Z+. Theorem 14.11 shows {z,a) = (0,0) is a 1:1 resonance point 
of (1491) . The numerical experiment is carried for m = 100, that is e = 

From Theorem 13.21 we know that 1^ — {(ai,a2) : |ai + |a2 + h.o.t. = 0} is the local 
representation of the Hopf point branch, while l^o = {(Q!i,a2) '■ ffo!! + + h.o.t. = 0} is 
the local representation of the homoclinic branch of (|48p in parameter plane {ai,a2). They 
are, by neglecting the higher order terms, plotted in Figure [1] by solid lines with triangles and 
diamonds, respectively. 

From Theorem l4.6l we know that /| = {(ai, a2) : (| + 2e)ai + (| - 2e)a2 + h.o.t. = 0} is the 
local representation of the Neimark-Sacker point branch, while l^^ = {(ai, a2) ■ (ff + ^£)oii + 
(ij — ^e)a2 + h.o.t. ~ 0} is the local representation of the homoclinic branch of the forward 
Euler method (j49p in parameter plane (ai, a2). They are, by neglecting the higher order terms, 
plotted in Figure [T] by solid lines with circles and stars, respectively. 

Besides, the realistic Hopf point branch of the forward Euler scheme (|49| is obtained by 
detecting the occurrence of the eigenvalues with modulus 1 at the fixed point z = 2{a2 — ai). It 
is plotted in Figure [T] by dash-dot line. The critical values of parameter a for the occurrence 
of homoclinic orbits of the forward Euler scheme are obtained by a shooting technique. 
See the dotted hue in Figure [TJ 

In the parameter plane (ai, a2), the fixed point of z = 2 (0:2 — ai) is a focus when (ai, 02) 
belongs to the left region of Z|, when the parameter (ai,Q!2) moves right and crosses it 
turns into a central point and there will be periodic solutions bifurcating from this point. The 
periodicity will tend to infinity as the parameter a moves right and tends to l'^. At last, the 
periodic solution becomes the homoclinic solution when a arrives at l^. These are the repro- 
duction of the bifurcation structures for (l48l) near the Takens-Bogdanov point. These processes 
are shown in Figure [2] to |4] for fixed a2 — —0.05, —0.15, —0.25, the corresponding values of ai 
are 0.005,0.05,0.08 (focuses), 0.028,0.085,0.145 (periodic solutions) and 0.0308,0.0950,0.1631 
(homoclinic solutions), where the derivatives of z are approximated by the difference quotient, 
and then a phase portrait like ODE's is presented. 
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1.18 



Figure 1: Local bifurcation diagram in parameter plane (01,02): by neglecting the h.o.t., the 
homoclinic branch loo and Hopf point branch Ih of pS)) . the theoretical homoclinic branch 
llo and Neimark-Sacker point branch If^ of the forward Euler discretization as well as the 
homoclinic branch "numerical HB" and neimark-S acker point branch "numerical HpB" detected 
in Forward Euler discretization ([^^ are plotted. 




(a) a = (0.005,-0.05) (b) a = (0.05, -0.15) (c) a = (0.08, -0.25) 

Figure 2: 2 (02 — ai) are focuses of the Euler method when a belongs to the left region of 

Il- 




ia.) a = (0.028, -0.05) (b) a = (0.085, -0.15) (c) a = (0.145, -0.25) 



Figure 3: Periodic solutions bifurcated from 2(q2 — ai) of the Euler method when a is 

located in the region confined by Z| and Z^. 
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0.6 




-0.6 - 

-0.8 ' ' ' ' ' ' ' 

-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 



Figure 4: Homoclinic orbits of (j49|): a = (0.0308,-0.05) gives the small one, a = 
(0.0950, —0.15) leads to the middle one, and a = (0.1631, —0.25) results in the large one. 
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